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Abstract:  Preconditioning  techniques  that  are  used  to  alleviate  numerical  stiffness 
due  to  low  Mach  numbers  in  steady  flows  have  typically  not  performed  well  for 
unsteady  low  Mach  problems  because  the  preconditioning  scaling  requirements  for 
preserving  discrete  accuracy  in  time-accurate  flows  are  very  different  from  those 
for  steady  flows.  Specifically,  distinct  scalings  are  necessary  for  the  velocity  and 
pressure  fields  under  the  low-Mach,  high-Strouhal  conditions  characteristic  of 
acoustic  wave  problems.  In  this  article,  a  unified  flux  formulation  is  presented 
where  the  optimal  scaling  required  for  spatial  accuracy  is  maintained  over  a  broad 
range  of  flow  conditions.  Both  upwind  flux-difference  and  AUSM-type  schemes 
are  investigated  and,  in  both  cases,  the  judicious  use  of  “steady”  and  “unsteady” 
preconditioning  scalings  in  the  flux  formulation  is  shown  to  be  critical  for 
preserving  accuracy.  Low  Mach  number  vortex  propagation  and  acoustic  problems 
are  used  to  demonstrate  the  strengths  of  the  formulation.  These  studies  show  that 
the  AUSM  family  generally  performs  better  than  the  blended  flux-difference 
schemes  especially  in  terms  of  vortex  shape  preservation. 

Keywords:  Low  Mach  Number  Preconditioning,  Unsteady  Flow. 


1  Introduction 

Accurate  and  efficient  modeling  of  unsteady  low  Mach  number  flows  continues  to  be  a  challenging 
problem  since  it  requires  the  resolution  of  disparate  time  scales.  Unsteady  effects  may  arise  from  a 
combination  of  hydrodynamic  effects  in  which  pressure  fluctuations  are  generated  by  unsteady 
velocity  fluctuations,  and  acoustic  effects  where  pressure  fluctuations  correlate  with  density 
fluctuations  even  if  the  mean  Mach  number  is  very  low.  Examples  in  the  former  category  would 
include  vortex  propagation  problems  and  bluff  body  shedding,  while  examples  of  the  latter  include 
both  aeroacoustic  problems  as  well  as  acoustic  wave  propagation  in  internal  flows  among  others. 
Many  practical  applications  including  rotorcraft  flows,  jets  and  shear  layers  include  a  combination  of 
both  acoustic  and  hydrodynamic  effects.  Furthermore  these  effects  may  be  localized  with  the  flow 
characteristics  varying  in  the  domain.  Typically  the  near-field  may  exhibit  non-linear  coupling 
between  the  different  modes  of  unsteadiness,  while  the  far-field  may  show  a  separation  of  vortex 
propagation  and  acoustic  scales  more  typical  of  low  Mach  number  flows.  Therefore,  it  is  important 
that  an  algorithm  designed  for  unsteady  low  Mach  number  flows  function  efficiently  over  a  broad 
range  of  flow  conditions  and  accurately  resolve  the  inherent  stiffness  of  the  system. 

Accurate  unsteady  simulations  of  these  complex  low  Mach  number  flows  pose  difficulties  for  time¬ 
marching  algorithms.  Preconditioning  techniques  that  are  traditionally  used  to  alleviate  numerical 
stiffness  work  well  for  steady  state  simulations  but  have  problems  with  both  efficiency  and  accuracy 
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for  unsteady  computations.  One  crucial  factor  arises  from  the  conflicting  requirements  placed  on  the 
flux  procedure:  the  preconditioning  scaling  that  is  optimal  for  accuracy  of  momentum  and  energy 
transport  is  different  from  that  needed  for  accuracy  of  the  acoustic  (pressure)  waves  under  the  low- 
Mach,  high-Strouhal  conditions  characteristic  of  acoustic  problems.  The  goal  of  the  current  effort  is 
the  development  of  a  unified  numerical  flux  procedure  where  the  optimal  scalings  required  for  spatial 
accuracy  of  all  fields  are  preserved  over  a  wide  range  of  flow  conditions.  A  related  aspect  of  this  work 
is  the  design  of  preconditioning  procedures  to  preserve  optimal  convergence  efficiency  as  well 
without  impacting  the  underlying  accuracy  of  the  discrete  formulation. 

The  development  of  preconditioning  algorithms  has  been  a  very  active  area  of  research  in  the 
literature  and  many  flavors  of  flux  procedures  have  been  developed.  These  schemes  can  be  broadly 
classified  into  two  families  based  on  their  general  characteristics  for  low  Mach  number  flows:  1) 
flux-difference/matrix-dissipation  preconditioning  algorithms  [l]-[4]  and  2)  the  AUSM  family  of 
flux-splitting  schemes  [5]-[8].  Both  sets  of  algorithms  have  been  applied  with  success  to  steady  low 
Mach  number  flows  by  recognizing  that  the  standard  flux  formulations  for  higher  Mach  numbers 
generate  too  little  dissipation  in  the  pressure  field  and  too  much  dissipation  in  the  velocity  field.  By 
reformulating  the  equations  by  introducing  a  preconditioning  matrix  (as  in  the  flux-difference 
schemes)  or  by  modifying  the  flux  formulation  using  a  preconditioning  scaling  term  (as  in  the 
AUSM+up  schemes),  the  stiff  acoustic  speeds  are  scaled  to  the  local  convective  velocity  and  good 
solution  convergence  and  accuracy  is  obtained  for  “all  speeds”. 

For  unsteady  flows  at  low  Mach  numbers  and  high  Strouhal  numbers  (needed  to  resolve  acoustic 
waves),  the  steady  preconditioning  scaling  that  effectively  filters  the  acoustics  from  the  solution 
becomes  far  too  dissipative  for  the  acoustic/pressure  field.  More  generalized  definitions  of  the 
preconditioning  parameter  have  been  proposed  that  take  the  local  Strouhal  number  into  account  [9]. 
Here,  the  preconditioning  parameter  reverts  back  to  its  non-preconditioned  value  as  acoustic  effects 
become  dominant  thus  restoring  the  ability  to  accurately  model  acoustic  propagation  with  good 
convergence.  However,  a  key  drawback  of  this  unsteady  preconditioning  scaling,  especially  for  the 
flux-difference  procedures,  is  that  the  formulation  becomes  inaccurate  for  the  velocity  and 
temperature  fields  which  are  still  governed  by  the  convective  fluid  time  scales  and  not  the  acoustic 
time  scales.  Therefore,  solutions  to  unsteady  low  Mach  number  problems  may  show  a  disconcerting 
disconnect  between  the  accuracy  requirements  of  the  velocity  and  the  acoustic  fields  [1]. 

In  an  effort  to  improve  on  this  behavior,  Sankaran  and  Merkle  reformulated  the  flux-difference 
scheme  with  features  akin  to  the  AUSM  schemes  [9].  Here,  the  eigenvalue  matrix  used  for  spatial 
flux  dissipation  was  altered  such  that  the  eigenvalue  associated  with  pressure  equation  corresponded 
to  the  acoustic  wave  from  the  unsteady  preconditioned  system  whereas  all  the  other  equations  use  the 
convective  velocity.  While  accurate  results  with  unsteady  preconditioning  are  reported  for  both 
acoustic  and  vortex  propagation  problem,  the  concern  with  this  formulation  is  that  the  spatial 
dissipation  does  not  reduce  to  the  “steady”  preconditioning  form  if  unsteady  effects  are  not  dominant. 
To  avoid  this  discrepancy  for  steady  calculations,  Potsdam  et  al.  [1]  formulated  a  “blended”  procedure 
by  defining  selection  matrices  that  use  the  flux  from  the  “unsteady”  preconditioning  matrix  for  the 
continuity  equation  (or  the  pressure  field)  and  the  flux  from  “steady”  preconditioning  for  the 
momentum  and  energy  equations  (or  the  velocity  and  temperature  fields)  [1].  In  the  limit  of  steady 
flow,  the  “unsteady”  preconditioning  parameter  for  the  continuity  equation  automatically  reverts  back 
to  the  steady  preconditioning  parameter  and  the  overall  steady  preconditioning  formulation  is 
recovered. 

The  second  class  of  numerical  flux  procedures,  the  AUSM  family  of  flux-split  schemes  has  also  been 
extended  to  low  Mach  number  flows  [5]-[8].  A  key  characteristic  of  AUSM  schemes  is  that  the 
convective  flux  is  separated  from  the  pressure  flux;  both  flux  terms  are  computed  independently  as 
scalar  formulations  thus  making  it  possible  to  independently  tailor  the  dissipation  for  hydrodynamic 
and  acoustic  unsteadiness.  The  convective  flux  terms  have  essentially  a  upwinded  form  with 
additional  dissipation  at  low  Mach  numbers;  Liou  demonstrates  excellent  results  for  steady  low  Mach 
number  flows  in  the  AUSM+up  formulation,  wherein  the  “u”  and  “p”  refer  to  velocity  and  pressure 
dissipation  terms  that  have  been  introduced  into  the  basic  AUSM+.  Specifically,  the  magnitude  of  the 
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pressure  dissipation  term  is  controlled  by  the  introduction  of  the  steady  preconditioning  scaling  in 
order  to  obtain  uniform  accuracy  of  this  term  at  all  Mach  numbers  [5].  Shima  and  Kitamura  [7], [8] 
have  proposed  a  variant  to  the  AUSM  scheme,  referred  to  as  the  SLAU  scheme  (Simple  Low 
Dissipation  AUSM),  which  has  been  applied  primarily  to  unsteady  flows.  Interestingly,  this  scheme 
does  not  involve  any  explicit  modifications  for  low  Mach  numbers,  but  as  we  will  later  show,  this 
scheme  is  naturally  formulated  for  unsteady  low-Mach  acoustic  problems. 

The  focus  of  this  paper  is  to  present  a  more  generalized  preconditioning  framework  based  on  an 
unsteady  Mach  number  parameter  which  ensures  that  the  flux  formulation  is  accurate  and  efficient  for 
both  hydrodynamic  and  acoustic  instabilities,  and  reverts  to  the  traditional  steady  flux  form  in  the 
limit  of  steady  low  Mach  number  flows.  The  generalized  unsteady  preconditioning  framework  has 
been  adapted  for  both  AUSM  and  flux-difference  schemes.  The  flux-difference  procedures  are 
derived  by  developing  “blended”  formulations  of  the  “unsteady”  and  “steady”  preconditioning 
parameters  and  are  extensions  of  the  work  presented  by  Potsdam  et  al.  [1].  For  the  AUSM  family  of 
schemes,  the  dissipation  parameters  for  both  the  mass  flux  and  pressure  terms  have  been  modified 
using  the  unsteady  Mach  number  parameter.  The  modified  AUSM  schemes  are  referred  to  as 
AUSM+up’  and  AUSM+u’p’,  where  the  “primes”  designate  that  the  pressure  and/or  velocity 
dissipation  terms  are  being  scaled  by  the  unsteady  Mach  number.  Corresponding  modifications  are 
proposed  for  the  SLAU  scheme  as  well.  Both  the  blended  flux-difference  and  the  modified 
AUSM+up  formulations  are  tested  on  a  range  of  unit  problems  encompassing  both  hydrodynamic  and 
acoustic  instabilities.  Solution  accuracy  and  unsteady  inner  iteration  convergence  are  evaluated  for 
these  problems  to  demonstrate  their  wide  range  of  applicability. 


2  Flux  Formulations  for  Low  Mach  Number  Unsteady  Flows 

The  formulation  of  a  more  generalized  unsteady  preconditioning  framework  for  both  flux  differenced 
and  AUSM  family  of  schemes  is  described  here.  We  begin  by  giving  a  brief  overview  of  the 
preconditioned  system  of  equations  and  steady  preconditioning  for-flux  difference  schemes  in 
Sections  2.1  and  2.2.  The  “blended”  flux  difference  formulations  for  unsteady  preconditioning  that 
were  developed  as  part  of  this  effort  are  described  in  Section  2.3.  The  AUSM  family  of  schemes  is 
described  in  Section  II. D  and  an  analysis  comparing  AUSM+up  and  SLAU  is  discussed  here.  The 
extensions  for  AUSM+up’  and  AUSM+u’p’  employing  the  unsteady  preconditioning  parameter  are 
described  in  Section  2.5. 

2.1  Preconditioned  Equation  System 

The  standard  conservative  form  of  the  equations  (with  a  two-equation  turbulence  model)  may  be 
written  as 


dQ  dE  dF  dG  0  ^ 
—  + - + - + - =S+D, 

dt  dx  dy  dz 


(1) 


where  Q  =  [p,  pu,  pv,  pw,  e,  pk,ps]T  represents  the  vector  of  conserved  variables;  E ,  F  and  G  are  the 
flux  vectors;  S  represents  the  source  terms  (if  any);  and  Dv  represents  the  viscous  fluxes.  For  low 
Mach  number  flows  as  density  approaches  a  constant  the  conservative  variables  become  ineffective 
for  temporal  integration.  In  this  flow  regime,  the  primitive  variables  vector,  Qv  =  [p,  u,  v,  w,  T,  k,  s]1 , 
constitutes  a  particularly  effective  choice.  Primitive  variables  replace  the  density  by  pressure,  thereby 
avoiding  round-off  errors,  and  the  total  energy  by  temperature  which  is  more  compatible  with  the 
thermal  diffusion  terms  and  for  modeling  multi-species  flows  (as,  for  example,  the  combusted  exhaust 
plume  from  an  aircraft  engine). 


One  effective  way  for  expressing  a  general  iterative  method  is  through  a  dual-time  formulation.  Upon 
appending  a  pseudo-time  derivative  and  a  preconditioning  matrix,  Tp,  Eqn.  (1)  takes  the  form; 


^  dQv  dQ  dE  dF  dG  __ 

p  +  + + - + =  S  +  Dr 

p  dr  dt  dx  dy  dz 


(2) 
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The  numerical  characteristics  of  the  pseudo-time  in  this  equation  are  determined  by  the  eigenvalues  of 
the  matrix,  (dE  /  Sgv)J  which  are: 


A  =  diag 


U(  1  +  R)±  ft2  (1  -  R)1  +  4c  '2  (42  +  ?  +  l\ ) 


tU9U9U9U9U 


(3) 


where  R  represents  the  ration  of  the  square  of  the  pseudo-acoustic  speed  and  the  physical  acoustic 
speed  which  can  be  expressed  as  the  square  of  the  ratio  of  the  physical  Mach  number  to  a  pseudo- 
Mach  number  as 


R  {/¥>&+ PrQ-phpj)  |VY  (m\ 

{pp'PK + pT<\- Phi)  vcj  W) 

This  system  of  equations  can  be  made  well-conditioned  by  choosing  p'p  to  be  of  order 
p'p  «[<9(l/w2)],  thereby  resulting  in  a  pseudo-sound  speed,  cr^[0{u)\.  The  preconditioning 
parameter,  p'p ,  can  be  expressed  in  terms  of  the  ratio  of  pseudo-Mach  number  to  physical-Mach 
number  and  the  ratio  of  specific  heats: 


Pp 

Pp 


r 


f  M' 


\Mp ) 


+  lzl 

r 


(5) 


where  M2p  -  min  [max(M,2 ,  M] , M2in ),  l]  (6) 

Here,  Mt  is  the  local  Mach  number,  Mmin  is  a  user  specified  cut-off  Mach  number  to  preclude 
difficulties  at  stagnation  points,  and  Mu  is  an  “unsteady”  preconditioning  scaling.  This  unsteady 
preconditioning  term  can  be  related  to  the  global  Strouhal  number  for  the  problem  as: 


Mu 


L 

7TcAt 


(7) 


where  A  is  a  global  length  scale  and  At  is  the  time-step.  In  the  following  discussions,  Eqn.  (6)  without 
the  unsteady  preconditioning  scaling  is  referred  to  as  steady  preconditioning  formulation,  while  the 
full  version  of  Eqn.  (6)  is  called  the  unsteady  preconditioning  formulation. 


2.2  Flux-Difference  Algorithm:  Roe’s  Scheme  with  Steady  or  Unsteady 
Preconditioning 

The  inviscid  flux  at  cell  interfaces  is  calculated  with  an  approximate  Riemann  solver  given  the 
reconstructed  left  and  right  solution  states.  The  flux  difference/matrix-dissipation  procedure  based  on 
Roe’s  scheme  [11]  can  be  given  as: 

Fm  =  ~  |  F(Qv  ’  5)  +  F(Qy  >  h)  +  AF~(QV ,  n)  +  AF+  (Qv ,  n)}  (8) 

where  h  is  the  area-directed  normal  of  either  the  cell  face  or  the  dual  face  crossing  the  edge,  Q;  and 
are  the  left  and  right  primitive  state  variables,  and  Qv  denotes  the  Roe-averaged  variables.  The 
flux  differences  AF-  and  AF+  are  given  by: 

AF+(Qv,noi)  =  ^(fpRp(Ap+  |  Ap  |)Lp)(Q+  -Q")  (9) 

AF  (Qv,noi)  =  -t(f  pRp(Ap- 1  Ap  |)Lp)(Qy  -Q~)  (10) 

Here  rp,  Rp,  Lp,  and  Ap  are  the  preconditioning  matrix,  right  eigenvector  matrix,  left  eigenvector 
matrix,  and  diagonal  eigenvalue  matrix  Qvfor  the  preconditioned  system  computed  with  the  Roe- 
averaged  variables. 
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2.3  Flux-Difference  Algorithms:  Blended  Unsteady/Steady  Preconditioned 
Schemes 

The  standard  flux  difference  formulation  shown  above  (Eqns.  (9)  and  (10))  indicate  that  the  spatial 
dissipation  is  intimately  tied  to  the  time  scaling  defined  by  the  preconditioning  matrix.  Table  1  gives 
the  behavior  of  the  pressure  and  velocity  dissipation  terms  for  both  steady  low  Mach  numbers  and  for 
unsteady  low-Mach,  high-Strouhal  numbers.  In  the  steady  low  Mach  case,  the  use  of  no 
preconditioning  clearly  leads  to  ill-behaved  artificial  dissipation  terms  with  the  pressure  dissipation 
becoming  vanishingly  small  and  the  velocity  dissipation  becoming  unboundedly  large.  The  use  of  the 
steady  or  inviscid  preconditioning  in  Eqn.  (6)  clearly  alleviates  this  situation  and  both  the  pressure 
and  velocity  dissipation  terms  become  the  same  magnitude  as  the  dominant  physical  terms  in  the 
equations  of  motion.  The  situation  for  unsteady  flows  in  the  low-Mach,  high-Str  limit  is  more 
challenging.  It  can  be  observed  that  using  just  the  inviscid  preconditioning  in  the  flux  formulation 
leads  to  increased  dissipative  errors  in  the  pressure  field,  while  the  velocity  dissipation  is  well- 
behaved.  In  contrast,  when  the  unsteady  preconditioning  scaling  given  in  Eqn.  (6)  is  used,  the 
pressure  dissipation  is  well  behaved,  while  the  velocity  dissipation  becomes  too  large.  Thus,  no 
scaling  of  the  standard  flux-difference  formulation  naturally  preserves  the  discretization  accuracy  in 
the  unsteady  case. 


Table  1:  Normalized  scaling  of  the  pressure  and  velocity  dissipation  terms  in  the  steady  low- 
Mach  limit  and  the  unsteady  low-Mach,  high-Str  limit  for  different  preconditioning  scalings  in 

the  flux-difference  scheme. 


Formulation 

Steady  Low  Mach  Limit 

Unsteady  Low-' 

Mach  High-Str 

Pressure 

Velocity 

Pressure 

Velocity 

No  preconditioning 

O(M) 

0(1 /M) 

0(1) 

0(1 /M) 

Steady  Preconditioning 

0(1) 

0(1) 

0(1 /M) 

0(1) 

Uns.  Preconditioning 

0(1) 

0(1) 

0(1) 

0(1 /M) 

Blended  scheme 

0(1) 

0(1) 

0(1) 

0(1) 

Also  shown  in  Table  1  is  the  so-called  blended  formulation  which  is  constructed  so  that  “steady”  and 
“unsteady”  preconditioning  forms  can  be  used  for  specific  terms  in  the  momentum  and  energy 
equations,  thereby  controlling  the  pressure  and  velocity  dissipation  scaling  independently.  These 
formulations  are  devised  by  exploiting  the  form  of  the  algebraic  expressions  for  the  spatial  dissipation 
terms.  The  combination  of  the  spatial  dissipation  terms  in  Eqns.  (9)  and  (10)  may  be  written  as  a 
matrix,  Cp,  multiplied  by  the  change  in  the  primary  dependent  variable  vector: 

AF+(Qv,noi)  + AF-(Qv,noi)  =  Cp  (Q+  -Q~)  (1 1) 

where  the  combined  dissipation  matrix,  Cp,  can  be  written  out  as 


Cn 

C12 

C13 

C14 

C15 

Ci6 

C17 

iiCn  +C21 

uC12  +C22 

uCl3  +C23 

llCj4  +  C24 

uC15  +C25 

uC16 

uC17 

vCll  +C31 

wcll  +C41 

HCn  +C5i 

HC12  +C52 

HCi3+C53 

HC14  +  C54 

HCi5+C55 

hc16 

hc17 

kCn 

kC12 

kC!3 

kC[4 

kC15 

kC16 

kC17 

eCn 

J 

The  form  of  this  matrix  shows  that  the  entries  for  the  scalar  transport  equations,  such  as  those  for  the 
turbulence  scalars  or  species  equations  for  combusting  flows,  are  multiples  of  the  first  row  (the 
continuity  equation)  scaled  by  the  convected  property  of  that  equation.  Similarly,  the  entries  for  the 
momentum  and  energy  equations  contained  one  term  that  is  a  product  of  the  first  row  multiplied  by 
the  pertinent  local  scalar  (e.g.,  u ,  H ),  however,  they  also  include  an  additive  correction  term  which  is  a 
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function  of  the  eigenvalues.  Therefore,  the  form  of  this  matrix  suggests  “blended”  schemes  in  which 
steady  and  unsteady  preconditioning  can  be  applied  to  different  parts  of  this  matrix. 

One  choice  of  blended  scheme  is  constructed  by  decomposing  the  dissipation  matrix  into  two  terms  as 
given  by 

AF+  +  AF-  =  (cOt+C^)(q+-Q-)  (13) 

where  the  first  term  is  given  by  the  tensor  product  of  two  vectors  given  by 
C  =  (Cn  Cl2  C13  C14  C15  C16  Cl7)  and  o  =  (l  u  v  w  H  k  s)  and  the  second  term  is  a 

correction  matrix  .  This  matrix  structure  suggests  a  sophisticated  blending  strategy  in  which 

“steady”  and  “unsteady”  preconditioning  forms  could  be  used  for  specific  terms  in  the  momentum  and 
energy  equations  which  may  be  written  as: 

AF+  +  AF_  =|cO  +Cp)(Q+-Q-)  (14) 

Here  the  COT  terms  use  “unsteady”  (green  notation)  while  the  terms  could  use  “steady” 

preconditioning  (red  notation).  This  blending  scheme  is  equivalent  to  having  the  first  row  of  the  Cp 
matrix  with  unsteady  preconditioning  while  the  other  rows  will  have  a  combination  of  steady  and 
unsteady  preconditioning  terms.  The  last  row  of  Table  1  shows  that  such  a  scaling  strategy  results  in 
well-behaved  pressure  and  velocity  terms  in  both  the  steady  and  unsteady  limits. 

Other  choices  of  blended  schemes  are  possible,  but  have  been  observed  to  yield  similar  scaling  results 
and  so  are  not  investigated  further  here.  Also,  we  note  that  this  blended  scheme  is  similar  (but  not 
identical)  to  the  scheme  proposed  by  Potsdam  et  al.  [1]. 


2.4  AUSM  Type  Algorithms  for  Preconditioned  Systems 

Flux-splitting  schemes  provide  an  alternative  approach  to  flux  difference  schemes  for  calculating  the 
inviscid  flux  at  cell  interfaces.  These  procedures  are  attractive  since  they  provide  a  natural  means  to 
decouple  the  dissipation  for  the  momentum  and  pressure  flux  thus  allowing  more  flexibility  to  tailor 
the  dissipation  for  acoustic  and  hydrodynamic  instabilities  in  unsteady  flows.  Prior  to  describing  the 
extensions  to  an  unsteady  preconditioning  framework  we  analyze  the  differences  between  the 
AUSM+up  scheme  [5]  developed  by  Liou  and  the  SLAU  scheme  derived  from  AUSM  by  Shima  and 
Kitamura  [7], [8]  specifically  for  low  Mach  number  flows.  These  schemes  are  described  below. 

The  numerical  flux  for  flux-split  schemes  is  given  by 

rh  +  \m\  ^  m  —  \m 

F=  - l—t0++ - - — 

2  2 

O  =  (\,  u,v,  w,hT)T  (16) 

The  unique  formulation  of  the  AUSM  family  schemes  is  determined  by  the  choice  of  mass  flux 
function,  m  ,  and  average  pressure,  p  .  The  mass  flux  in  AUSM+up  scheme  is  given  by 


-0  +  pn 


(15) 


m  =  cMV2 


ifMm>  0 
otherwise  ’ 


MV2=M+w(Ml)  +  Mw(Mr) 


+^max(l  -  (jM2 , 0)  Pr_J' 

fa  PC 


(17) 


Here  M(4)  are  fourth  order  polynomials  in  Mach  number.  The  last  term  is  a  dissipation  term 

formulated  for  low  Mach  number  flows.  The  fa  term  in  the  denominator  increases  the  dissipation  in 
inverse  proportion  to  the  Mach  number  and  is  equivalent  to  a  steady  preconditioning  parameter.  It  is 
defined  as: 
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/a=Mo(2-Mo)e[0,l],  M2=min(l,max(M2,M2in)),  M2  =(u2L +u2R)/2c2  (18) 

Making  a  low  Mach  number  assumption  and  dropping  higher  order  terms  in  Mach  number  Eqn.  (17) 
can  be  rewritten  as: 


m  =  P±  I  -( Vn+  +  v; )  -  max(l  -aM2, 0)  Pr_?l  \  ( 1 9) 

1+  fa  pc  ) 

where  p±  is  p+  for  Mm  greater  than  zero  and  p~  otherwise. 

The  pressure  flux  formulation  in  the  AUSM+up  scheme  is  given  by; 

P  =  I&(ML)pL+ir(MR)pR -2KuP5+P5pfacAu  (20) 

Here  i**  are  fifth-order  polynomials  in  Mach  number  while  the  last  term,  which  can  be  summarized 

as  the  pu  dissipation  term,  is  a  dissipation  term  that  operates  at  transonic  Mach  numbers  and  decreases 
in  value  as  the  Mach  number  drops.  Making  a  low  Mach  number  assumption  and  dropping  higher 
order  terms  in  Mach  number  this  equation  reduces  to: 

p  --^iyP  +  p~^  +  ^(M+p+  -M~p~^-^[2M~\pcAu  (21) 


where  the  [2M]  factor  in  the  last  term  comes  from  the  expression  forfa. 

For  the  SLAU  scheme  given  by  Shima  and  Kitamura  [7], [8],  the  mass  flux  is  given  directly  in  the 
simplified  low  Mach  number  form  as 

f  =  -2{piv;+\V,[)  +  p-(v,--\v,\-)}-JL*p  (22) 

*  =  1-2  M 

While  the  pressure  flux  for  the  SLAU  formulation  is  given  by 


~pJ-£±£2+  V-jO+(l-;ri(% +%  (23) 

Making  a  low  Mach  number  assumption  and  dropping  higher  order  terms  in  Mach  number  this 
equation  reduces  to: 


P~)  +  \{M+  +  Ar)(+  ~P~)-\[2M]pc/±u 

O  O 


(24) 


where  the  [2M]  factor  in  the  last  term  comes  from  the  expression  for  %. 

Comparing  the  mass  flux  dissipation  for  AUSM+up  (Eqn.  (19))  and  SLAU  (Eqn.  (22)),  it  is  observed 
that  the  primary  difference  for  the  two  schemes  is  the  dissipation  that  is  added  at  low  Mach  numbers 
for  the  AUSM+up  scheme.  In  the  AUSM+up  scheme,  the  dissipation  term  for  the  mass  flux  (which  is 
a  pressure  dissipation)  has  the  fa  term  in  the  denominator  which  increases  the  dissipation  substantially 
as  the  Mach  number  drops  and  behaves  like  a  steady  preconditioning  parameter.  In  contrast,  the 
SLAU  formulation  does  not  have  this  term  and  exhibits  only  a  weak  dependence  on  Mach  number 
from  the  %  term  in  the  numerator.  Therefore,  the  SLAU  scheme  would  be  inadequate  for  steady  low 
Mach  number  problems  as  seen  in  Table  2,  which  shows  the  behavior  of  the  pressure  and  velocity 
dissipation  terms  for  the  steady  low-Mach  and  unsteady  low-Mach/high-Str  limits. 

In  contrast,  the  AUSM+up  scheme  is  well-behaved  for  steady  low-Mach  problems.  Interestingly,  in 
the  unsteady  acoustic  limit,  the  SLAU  scheme  (which  is  similar  to  the  basic  AUSM+  scheme) 
performs  well,  while  the  AUSM+up  scheme  introduces  too  much  damping  in  the  acoustic/pressure 
field.  To  obtain  uniform  accuracy  under  both  steady  and  unsteady  limits,  we  consider  in  the  following 
section  the  extensions  of  the  unsteady  preconditioning  to  the  AUSM  family  of  schemes. 
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Table  2:  Normalized  scaling  of  the  pressure  and  velocity  dissipation  terms  in  the  steady  low- 
Mach  limit  and  the  unsteady  low-Mach,  high-Str  limit  for  different  preconditioning  scalings  in 

the  AUSM  family  of  schemes. 


Formulation 

Steady  Low  Mach  Limit 

Unsteady  Low-Mach  High-Str  Limit 

Pressure 

Velocity 

Pressure 

Velocity 

SLAU 

O(M) 

0(1) 

0(1) 

0(1) 

AUSM+up 

0(1) 

0(1) 

0(1 /M) 

0(1) 

AUSM+up’ 

0(1) 

0(1) 

0(1) 

0(1) 

AUSM+u’p’ 

0(1) 

0(1) 

0(1) 

0(1 /M) 

2.5  Extensions  of  Unsteady  Preconditioning  Framework  to  AUSM  schemes 

As  our  analysis  in  the  previous  section  indicates  the  standard  AUSM +up  scheme  is  optimized  for 
steady  low  Mach  number  problems.  A  more  general  unified  formulation  that  automatically  selects  the 
appropriate  dissipation  level  is  presented  here  by  formulating  the  mass  flux  dissipation  in  terms  of  an 
unsteady  Mach  number  scale.  This  modified  scheme  is  referred  to  AUSM +up’\  where  the  pressure 
prime  indicates  that  the  pressure  dissipation  is  scaled  by  the  unsteady  Mach  number.  The  dissipation 
for  the  mass  flux  term  in  the  AUSM +up  ’  scheme  is  generalized  by  modifying  the  definition  of  Mach 
number  used  to  define  the  fa  term  in  Eqn.  (18)  by  including  an  unsteady  preconditioning  scale  as 
follows: 


M2o  =  min  (l,  max  (m1  ,  M^m  ,M„2))  (25) 

We  note  that  in  the  most  general  configuration  the  Mu  scale  can  be  defined  both  as  a  global  parameter 
and  as  local  factor  that  is  computed  on  the  local  flow  physics.  For  simpler  problems  a  global 
parameter  might  suffice  but  for  more  complex  physics  where  the  characteristics  may  vary  (e.g.  high 
frequency  in  the  near  field  of  a  turbulent  jet  and  low  frequency  in  the  far-field)  the  ability  for  the 
numerical  formulation  to  select  the  automatically  select  the  unsteady  scale  is  crucial.  We  note  that 
reducing  the  dissipation  for  the  mass  flux  term  for  unsteady  acoustic  flows  has  been  suggested  both 
by  Liou  [5]  and  Vigneron  et  al.  [6]  however  Eqn.  (25)  provides  a  means  for  the  numerical  formulation 
to  self-select  the  appropriate  dissipation  level.  Table  2  confirms  that  the  resulting  scheme  is  indeed 
well-behaved  in  both  the  steady  and  unsteady  limits  of  interest. 

An  additional  modification  involves  changing  the  pressure  dissipation  term  by  adding  additional 
dissipation  for  low  Mach  number  acoustic  problems  using  the  unsteady  Mach  number  parameter. 
These  modifications  lead  to  the  AUSM+w  ’p  ’  scheme  where  the  primes  indicate  that  both  pressure  and 
velocity  dissipation  terms  are  scaled  by  the  Mach  number  parameter.  The  definition  of  Mach  number 
M0  using  the  unsteady  Mach  number  scale  Mu  as  shown  above  in  Eqn.  (25)  also  affects  the  pressure 
dissipation  term  in  Eqn.  (20)  as  shown  below: 

~  0  when  AT  «  0  for  steady  low  Mach  number  flows 
a  u  (26) 

fa  « 1  when  Mu  « 1  for  unsteady  low  Mach  number  flows 

In  particular,  for  steady  flows  no  additional  dissipation  results  but  for  unsteady  acoustic  low  Mach 
number  flows  substantial  dissipation  is  added  to  the  pressure  flux.  Table  2  shows  that  such  a 
formulation  adversely  impacts  solution  accuracy  in  the  velocity  field  for  unsteady  acoustic  problems. 
Thus,  the  AUSM+up’  formulation  appears  to  be  the  ideal  choice  of  scheme  for  a  wide  range  of  steady 
and  unsteady  flow  conditions  with  the  exception  of  a  narrow  class  of  acoustic  propagation  problems 
where  the  AUSM+w'/?'  was  required  (see  discussion  in  Section  3.4).  We  further  note  that  the 
dissipation  flux  modifications  in  AUSM+up,  AUSM+up’  and  AUSM+u’p’  can  be  implemented  in  the 
context  of  the  SLAU  scheme  as  well  with  similar  scaling  behaviors  of  the  associated  dissipation 
terms. 
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3 Evaluation  of  the  Low  Mach  Number  Formulations 


The  low  Mach  number  formulations  described  in  the  previous  section  were  evaluated  with  several  test 
problems  involving  hydrodynamic  and  acoustic  unsteadiness.  These  calculations  were  performed 
with  the  CRUNCH  CFD®  code,  developed  at  CRAFT  Tech  [12]-[15],  The  candidate  flux 
formulations  for  unsteady  low  Mach  number  flows  will  be  tested  out  rigorously  for  the  following 
three  test  cases  that  include  both  hydrodynamic  and  acoustic  instabilities:  1)  Unsteady  inviscid  Lamb 
vortex  problem  (hydrodynamic  instability),  2)  Unsteady  inviscid  flow  in  a  pipe  with  fluctuating  back 
pressure  (mixed  acoustic  and  hydrodynamic  instability),  and  3)  Shock  tube  with  small  pressure 
difference  (pure  acoustic  problem).  We  note  that  both  the  “blended”  flux  difference  and  AUSM  type 
schemes  have  to  be  run  with  inconsistent  LHS  and  RHS  discretization  since  these  formulations  are 
unstable  when  a  consistent  LHS  is  used.  For  the  results  presented  here  we  employ  the  standard 
preconditioned  Roe  flux  differenced  procedure  on  the  LHS  with  unsteady  preconditioning  and  dual 
time  stepping  at  each  time  step.  Prior  to  discussing  the  unsteady  test  cases,  a  steady  low  Mach  number 
test  case  is  presented  to  compare  the  Roe  flux  differenced  procedure  with  the  SLAU  scheme;  as 
discussed  earlier  the  dissipation  in  the  original  SLAU  scheme  is  shown  to  be  appropriate  only  for 
unsteady  flows  and  this  constraint  is  demonstrated  by  computing  a  steady  flow-field. 


(b)  Pressure  coefficient  (c)  Residual  convergence 

Figure  1.  Solution  for  a  steady  flow  over  a  NACA0015  airfoil  at  a  freestream  Mach  number  of 
0.001  and  an  angle  of  attach  of  4°  as  computed  with  the  SLAU  scheme  and  a  Roe  flux 
differenced  scheme  with  steady  preconditioning. 
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3.1  Steady  Flow  over  a  NACA0015  Airfoil 

Steady  state  calculations  were  performed  for  a  NACA0015  airfoil,  at  a  freestream  Mach  number  of 
0.001,  Reynolds  number  of  1.95E+06,  and  an  angle-of-attack  of  4°  using  a  Roe  flux  differences 
procedure  with  steady  preconditioning  and  the  SLAU  scheme.  The  higher  order  solution  and 
convergence  is  shown  Figure  1  for  both  schemes.  For  the  SFAU  scheme  it  is  noted  that  the  solution 
shows  oscillations  with  odd-even  coupling  that  is  reflected  in  the  surface  pressure  coefficient. 
Moreover,  the  corresponding  convergence  shown  in  Figure  1(c)  indicates  that  the  SFAU  solution  does 
not  converge  in  contrast  to  the  solution  with  flux  difference  procedure  with  steady  preconditioning. 
The  numerical  instability  for  the  SFAU  was  verified  to  result  from  the  low  freestream  Mach  numbers; 
the  odd-even  oscillations  were  not  evident  when  the  Mach  number  was  increased  to  0.1  (results  not 
shown).  Therefore,  it  was  concluded  that  the  mass  flux  dissipation  in  Eqn.  (22)  for  SFAU  scheme 
was  unable  to  provide  accurate  solutions  for  steady  low  Mach  number  flows  as  expected. 

The  dissipation  term  for  SFAU  was  subsequently  modified  by  adding  the  form  of  the  dissipation  from 
the  AUSM+up  formulation  (Eqn.  (19)).  Therefore,  for  steady  low  Mach  number  problems  substantial 
dissipation  (resulting  from  the  fa  term  in  the  denominator  of  Eqn.  (19)  is  added  to  the  mass  flux  term, 
leading  to  what  we  refer  to  as  the  SFAU+p  scheme.  The  solution  for  the  NACA0015  airfoil  with  the 
modified  dissipation  term  is  shown  in  Figure  2.  The  flow  contours  and  pressure  coefficient  are  now 
smooth  and  the  convergence  was  found  to  be  independent  of  the  Mach  number.  Therefore,  this 
confirms  our  premise  discussed  earlier  that  the  SFAU  scheme  corresponds  to  the  standard  AUSM+ 
formulation  and  requires  modifications  in  the  steady  limit  (i.e.,  the  SFAU+p  scheme)  for  low  Mach 
number  flows. 


(a)  Pressure  contours  with  additional  dissipation  at  low  Mach  number  number  (SLAU+p  scheme) 


(b)  Pressure  coefficient  (c)  Residual  convergence 

Figure  2.  Solution  for  a  steady  flow  over  a  NACA0015  airfoil  at  a  freestream  Mach  number  of 
0.001  and  an  angle  of  attach  of  4°  as  computed  with  the  SLAU+p  scheme. 
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3.2  Vortex  Propagation  Problem 

The  propagation  of  an  unsteady  convecting  inviscid  Lamb  vortex  was  used  to  assess  the  effectiveness 
of  the  flux  formulations  discussed  above  for  the  hydrodynamic  class  of  problems.  The  velocity 
distribution  in  polar  coordinates  for  the  Lamb  vortex  is  given  by 


Vr  =  0 


v0=r 


l-e 


J.  ,j2 


(27) 


v  j 

where  r and  (j)  are  the  vortex  strength  and  characteristic  radius  of  the  vortex  and  were  set  to  O.lMoo 
and  0.03,  respectively.  Simulations  are  presented  here  for  a  freestream  Mach  number  of  0.001  were 
computed  on  a  uniform  Cartesian  grid  with  a  spacing  of  0.005.  The  vortex  solution  is  presented  after 
the  vortex  has  travelled  a  distance  equal  to  0.4  on  this  uniform  mesh.  The  unsteady  Mach  number 
scale  was  specified  as  done  by  Potsdam  et  al.  in  Ref.  [1]. 


M  =  L(MJ  (28) 

u  tt(Ax)(CFLu ) 

For  a  baseline  value  of  CFLU  =  1  this  translates  to  an  unsteady  Mach  number  of  2.55E-02.  For  very 
small  time  steps  (e.g.,  CFLU  =  0.001)  the  unsteady  Mach  number  becomesl  and  the  dissipation  levels 
are  representative  of  unsteady  preconditioning.  Results  are  presented  for  the  baseline  time-step 
dictated  by  CFLU  =  1  and  an  extremely  small  time-step  given  by  CFLU  =  0.001.  Two-hundred  and 
fifty  inner  iterations  were  used  for  all  cases.  Both  the  blended  flux  difference  scheme  and  the 
AUSM+up’  scheme  are  compared  with  benchmark  values  of  Roe  flux  differenced  scheme  using 
steady  and  unsteady  preconditioning. 

The  vortex  solution  and  inner  iteration  convergence  are  shown  for  the  baseline  time  step  of  CFLU  =  1 
in  Figure  3  for  the  different  schemes.  The  exact  vorticity  profile  is  given  by  the  black  line  in  Figure 
3(a).  The  red  and  blue  lines  correspond  to  the  standard  flux-difference  scheme  with  steady  and 
unsteady  preconditioning,  respectively.  The  result  computed  with  a  blended  flux  formulation  is  given 
by  the  green  line  whereas  the  orange  line  gives  the  profile  computed  with  the  AUSM+up’  scheme. 
The  vorticity  profiles  show  that  the  flux  difference  scheme  with  steady  preconditioning,  the  blended 
flux  difference  scheme,  and  the  AUSM+up’  collapse  to  essentially  the  same  solution.  The  solution 
computed  using  the  flux  difference  scheme  with  unsteady  preconditioning  has  poor  accuracy 
compared  to  the  other  schemes,  however,  the  convergence  plots  shows  that  it  has  the  best  rate  of 
convergence  indicating  optimal  time  scaling  for  convergence.  Rectification  of  this  inconsistency 
between  solution  accuracy  and  convergence  was  one  of  the  primary  goals  at  the  outset  of  this 
research.  The  convergence  for  the  blended  scheme  is  slightly  better  than  the  flux  difference  scheme 
with  steady  preconditioning.  The  convergence  for  the  AUSM+up’  scheme  is  worse  than  the  blended 
scheme  which  is  likely  due  to  the  increased  inconsistency  of  the  left-hand  side  and  right-hand  side 
flux  formulations,  however,  the  solution  was  just  as  accurate  as  the  blended  and  pure  steady 
preconditioning  schemes. 


11 

Distribution  Statement  A:  Approved  for  public  release;  distribution  is  unlimited. 


(a)  Vorticity  profile 


(b)  Inner  iteration  convergence 


Figure  3.  Vortex  solution  accuracy  and  convergence  a  time-step  of  CFLU  =  1. 


To  further  confirm  the  insensitivity  of  the  solution  accuracy  to  the  unsteady  preconditioning 
parameter,  the  time-step  was  reduced  to  CFLU  =  0.001  which  goes  to  the  limit  of  unsteady 
preconditioning  for  the  blended  and  AUSM+up’  schemes.  The  vorticity  profiles  and  inner  iteration 
convergence  are  plotted  in  Figure  4.  At  this  much  smaller  time  step,  the  steady  preconditioning  case 
and  AUSM+up’  show  essentially  flat  convergence  while  the  baseline  unsteady  preconditioning  and 
blended  scheme  show  rapid  convergence.  From  Figure  4(a),  it  is  observed  that  the  accuracy  for  the 
baseline  unsteady  preconditioning  scheme  deteriorates  very  dramatically  while  the  blended  scheme 
continues  to  provide  an  accurate  vorticity  profile.  Despite  the  very  poor  convergence,  the  AUSM+up’ 
scheme  and  the  baseline  scheme  with  steady  preconditioning  give  as  accurate  a  profile  as  the  blended 
flux  formulation. 


(a)  Vorticity  profile 


(b)  Inner  iteration  convergence 


Figure  4.  Vortex  solution  accuracy  and  convergence  for  a  time-step  of  CFLU  =  0.001. 


Thus  far  the  predicted  vorticity  field  was  considered  as  the  metric  for  evaluating  the  solution 
accuracy.  However,  while  the  pressure  variations  may  be  very  small,  its  gradient  should  correspond 
to  the  swirl  velocity  in  the  vortex.  Therefore,  it  is  instructive  to  evaluate  the  resulting  pressure 
contours  for  the  difference  flux  procedures.  The  vorticity  and  pressure  contours  for  the  various 
schemes  are  plotted  in  Figure  5  and  Figure  6,  respectively.  The  flux  difference  scheme  with  steady 
preconditioning  gives  accurate  results  for  the  vorticity  field  as  seen  in  Figure  5(a),  however,  the 
pressure  field  in  Figure  6(a)  is  very  inaccurate.  This  is  not  unexpected  since  the  steady 
preconditioning  gives  the  correct  scale  for  the  velocity  equation  but  provides  an  incorrect  scale  for  the 
acoustic  eigenvalues.  The  flux  difference  scheme  with  unsteady  preconditioning  gives  poor  solutions 
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for  both  the  vorticity  and  pressure  fields  in  Figure  5(b)  and  Figure  6(b)  due  to  large  dissipation  in  the 
velocity  equation. 


(c)  Blended  Scheme 


(d)  AUSM+wp ' 
Scheme 


(a)  Flux  Difference  (b)  Flux  Difference 
with  Steady  with  Unsteady 

Preconditioning  Preconditioning 

Figure  5:  Vorticity  field  for  various  flux  schemes  for  a  time  step  of  CFLU=0.001. 


The  blended  flux  difference  scheme  gives  the  correct  solution  for  both  the  vorticity  and  pressure  field 
as  shown  in  Figure  5(c)  and  Figure  6(c),  since  they  use  different  preconditioning  scales  for  the  two 
variables:  unsteady  for  the  pressure  wave  and  steady  for  the  velocity  field.  However,  there  appears  to 
be  one  notable  discrepancy  in  that  the  pressure  field  is  not  circular  to  match  the  vortex  but  is  instead 
distorted  to  be  a  rhombus.  The  cause  of  this  distortion  is  not  fully  understood  but  may  be  due  to  the 
difference  in  magnitude  for  the  cross-dissipation  terms  involving  the  velocity  components  and  the 
pressure  field  (e.g.,  AuAp  versus  AvA p). 


The  vorticity  and  pressure  fields  computed  using  the  modified  AUSM+up’  scheme  can  be  seen  in 
Figure  5(d)  and  Figure  6(d).  Accurate  solutions  are  realized  for  both  the  velocity  and  pressure  fields. 
In  particular,  the  pressure  field  does  not  show  the  distortion  that  was  evident  from  the  blended  scheme 
results.  Moreover,  the  circular  contours  in  the  pressure  field  accurately  match  the  contours  of  the 
velocity  field.  The  blended  flux-difference  and  the  AUSM+up’  schemes  both  show  odd-even 
oscillations  in  the  pressure  field.  The  source  of  this  instability  is  not  clear  but  may  be  related  to  the 
inconsistency  between  the  schemes  used  for  the  flux  for  the  implicit  and  explicit  sides. 


(c)  Blended  Scheme 


(d)  AUSM+up ' 
Scheme 


(a)  Flux  Difference  (b)  Flux  Difference 
with  Steady  with  Unsteady 

Preconditioning  Preconditioning 

Figure  6:  Gauge  pressure  field  for  various  flux  schemes  for  a  time  step  of  CFLU=.001. 


In  summary,  both  the  blended  flux  difference  and  the  AUSM+up’  schemes  are  able  to  provide 
accurate  solutions  of  the  vorticity  field  and  are  insensitive  to  the  time  step.  Both  schemes  are  a 
substantial  improvement  over  previous  baseline  flux  difference  procedures  that  showed  inconsistency 
between  convergence  and  accuracy  with  the  solution  accuracy  in  particular  deteriorating  rapidly  with 
unsteady  preconditioning.  An  evaluation  of  the  pressure  contours  revealed  that  both  the  blended  flux 
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difference  and  the  AUSM+up’  schemes  were  able  to  provide  the  correct  pressure  depression  in  the 
vortex.  However,  the  solutions  for  the  blended  flux  difference  scheme  show  a  distortion  of  the 
pressure  field  shape.  In  contrast,  the  pressure  field  computed  using  the  AUSM+up’  scheme 
maintained  the  circular  pressure  field.  The  only  drawback  of  the  AUSM+up’  formulation  appears  to 
be  less  than  ideal  convergence  at  small  time  steps  and  this  is  probably  related  to  the  inconsistent 
discretizations  on  the  LHS  and  RHS.  This  is  potentially  an  area  requiring  improvement  in  terms  of 
efficiency  of  the  procedure. 


3.3  Simulations  of  Flow  in  a  Tube  with  Oscillating  Back  Pressure 

The  next  test  problem  investigated  was  inviscid  flow  in  a  1-D  tube  where  the  back  pressure  was 
fluctuated  and  the  inlet  pressure  was  held  constant  as  given  by 


px=L=Pb(  1  +  ^sinwO 
P°=0  =  constant 

The  analytical  solution  for  this  scenario  is  given  below: 


(29) 


u\t)  =  - 


\Plu o  j 


vl  +  Q2 

P'(x,  t)  =  [s  sin  wt  +  j ou0u] 


sin  wt  -  Q  cos  wt  +  Qe  L  , 

—  -  y ou0u\  Q  =  , Phase  lag</>  =  tan”1  (Q) 

L  u0 


(30) 


The  velocity  field  is  uniform  in  space  but  fluctuates  in  time,  while  the  pressure  field  has  a  linear  slope 
spatially  from  the  inlet  to  the  exit  while  also  fluctuating  in  time.  The  key  parameter  of  interest  here  is 
the  Strouhal  number  Q;  as  Q  increases  the  unsteady  time  scales  become  dominant  and  the  baseline 
flux  difference  procedures  provide  inaccurate  solutions  that  vary  with  the  preconditioning  parameter. 
Therefore,  this  problem  would  be  a  good  test  for  the  validity  of  both  the  blended  flux  difference  and 
the  modified  AUSM Pup  ’  schemes.  The  objective  is  to  be  able  to  get  an  accurate  solution  independent 
of  the  preconditioning  parameter  and  frequency. 


In  the  present  calculations,  the  mean  Mach  number  in  the  tube  was  specified  to  be  0.001  and  a  5 
percent  fluctuation  amplitude  was  applied  at  the  exit  with  100  inner  iterations  for  the  dual  time 
iteration.  The  frequency  and  the  time  step  were  varied  over  a  broad  range  to  test  the  robustness  and 
accuracy  of  the  flux  procedures  over  a  broad  range  of  conditions.  The  steady  preconditioning  Mach 
number  cut-off  was  specified  as  0.001  while  the  unsteady  Mach  number  scale  was  set  to  1.0.  The 
time  step  was  specified  as  2  n/(PPW Q)  where  PPW denotes  the  points  per  time  period. 


For  the  benign  scenario  where  the  fluctuation  frequency  is  low,  Q  =  1,  and  the  time-step  is  large  with 
10  points  per  time  period,  the  baseline  flux  difference  scheme  with  steady  and  unsteady 
preconditioning,  the  blended  flux  difference  scheme,  and  the  AUSM+up’  scheme  all  produced 
pressure  and  velocity  profiles  that  matched  the  exact  solution.  Essentially  in  this  case  the  flow 
behaves  like  a  quasi-steady  problem  with  the  time  variation  being  a  series  of  steady  variations.  The 
convergence  history  shows  very  little  difference  between  steady  and  unsteady  preconditioning  for  the 
various  schemes  and  is  consistent  with  this  quasi-steady  picture. 


The  results  for  the  more  difficult  case  in  which  the  fluctuation  frequency  is  high,  Q  =  100,  and  the 
time-step  is  small  with  1000  points  per  time  period  are  presented  in  Figure  7.  For  the  flux  difference 
scheme  with  steady  preconditioning,  the  pressure  response  becomes  very  unstable  and  the  amplitude 
overshoots  the  exact  solution  for  velocity  and  pressure  as  seen  in  Figure  7(a)  and  Figure  7(b), 
respectively. 
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(c)  Velocity  residual  (d)  Pressure  residual 

Figure  7:  Velocity  and  pressure  history  and  residuals  for  a  high  oscillation  frequency 
(Q  =  100 )  computed  with  a  small  time  step  (1000  PPW) 


Moreover,  an  excessive  phase  error  is  also  evident.  The  pressure  residual  for  the  steady 
preconditioning  scheme  is  shown  in  Figure  7(d)  and  shows  no  convergence  during  the  inner 
iterations.  This  is  consistent  with  the  incorrect  flow  results.  The  solutions  computed  using  the 
blended  flux  difference  scheme,  and  the  AUSM+up’  gave  essentially  the  same  results  that  matched 
the  exact  solution  very  well.  Therefore,  only  the  solution  profiles  predicted  by  the  modified 
AUSM+up’  scheme  are  shown  in  Figure  7.  The  velocity  profile  matches  the  exact  solution  perfectly. 
However,  an  initial  fluctuation  around  the  exact  solution  can  be  seen  for  the  pressure  transient  in 
Figure  7(b).  This  fluctuation  decays  very  quickly  and  the  source  of  the  error  is  not  clear  at  this  point. 
The  convergence  history  plotted  in  Figure  7(c)  and  Figure  7(d)  for  the  velocity  and  pressure  residuals 
show  that  the  AUSM+up’  show  good  convergence  slopes  in  general  for  both  the  velocity  and  pressure 
residual.  In  contrast  the  steady  preconditioning  shows  very  poor  convergence  consistent  with  the  fact 
that  steady  preconditioning  is  inappropriate  for  high  frequency  pressure  fluctuations.  These  results 
confirm  that  for  low  Mach  number  acoustic  problems  that  the  dissipation  for  the  continuity  equation 
(pressure  variable)  for  both  the  AUSM+up’  and  the  blended  flux  difference  form  must  be  consistent 
with  the  unsteady  preconditioning  form. 

3.6  Simulations  of  Shock  Tube  with  Small  Pressure  Difference 

The  final  test  case  is  a  pure  acoustic  problem  with  low  velocities.  Here  we  model  a  shock  tube  with  a 
very  small  pressure  difference  across  the  diaphragm  that  generates  a  low  Mach  number  flow.  In  this 
particular  case,  a  tube  of  length  1  m  was  chosen  and  the  diaphragm  was  placed  at  x  =  0.5.  The 
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pressure  on  the  left  side  of  the  diaphragm  is  initialized  to  100028.04  Pa  and  the  pressure  on  the  right 
side  is  initialized  to  100000  Pa;  therefore  the  pressure  difference  is  28.04  Pa.  The  temperature  is 
initialized  to  300  K  on  both  sides  and  the  initial  velocity  is  zero.  A  very  weak  acoustic  waves  results 
and  generates  a  low  Mach  number  flow  of  the  order  of  0.0001  in  the  contact  interface. 

In  the  calculations  presented  here,  the  time  step  was  kept  constant  at  10  ps  and  the  calculations  were 
done  for  a  total  time  of  750  ps.  One -hundred  inner  iterations  were  used  at  each  time-step  and  the 
unsteady  Mach  number  for  unsteady  preconditioning  was  set  to  1.0.  The  results  of  the  previous 
problem  confirmed  that  dissipation  corresponding  to  unsteady  preconditioning  must  be  employed 
when  simulating  low  Mach  number  acoustic  problems.  Therefore,  the  flux  difference  procedure  with 
steady  preconditioning  will  not  be  considered  here. 

Figure  8  shows  the  solution  predicted  using  the  baseline  flux  difference  scheme  with  unsteady 
preconditioning  and  the  blended  flux  difference  procedure.  Recall  that  the  blended  scheme  uses 
unsteady  preconditioning  for  the  pressure  waves  and  steady  preconditioning  for  the  velocity  and 
temperature  equation.  For  the  solution  profiles  shown  in  Figure  8,  the  steady  preconditioning 
parameter  is  set  to  Ms  =  0.01  and  the  unsteady  Mach  number  scale  was  set  to  Mu  =  1 .0.  It  can  be  seen 
in  Figure  8(a)  and  Figure  8(b)  that  the  acoustic  wave  and  the  contact  interface  have  propagated  to  the 
correction  location  and  the  inner  iteration  convergence  history  shown  in  Figure  8(c)  appears  to  be 
good  for  both  flux  difference  procedures. 


Figure  8:  Shock  tube  solutions  for  the  blended  flux  difference  schemes  with  Mu=1.0  and 

Ms=0.01 


However,  closer  inspection  shows  that  the  blended  flux  difference  procedure  produces  small 
fluctuations  in  the  pressure  field  at  the  contact  interface  which  translates  to  substantial  fluctuations  in 
the  velocity.  The  level  of  the  fluctuations  is  sensitive  to  the  disparity  between  the  steady  and  unsteady 
preconditioning  in  the  blending  formula;  the  oscillation  amplitude  increases  as  the  steady  cut-off 
Mach  number  is  reduced  and  the  oscillation  amplitude  decreases  as  the  steady  cut-off  Mach  number  is 
raised.  If  Ms  is  increased  to  1.0  then  the  blended  flux  procedure  is  identical  to  the  flux  difference 
scheme  with  unsteady  preconditioning  which  does  not  produce  any  oscillations. 

The  profiles  predicted  using  the  AUSM +up  schemes  are  shown  in  Figure  9  for  two  different 
formulations.  The  profiles  given  by  the  blue-line  include  the  modification  to  the  mass  flux  dissipation 
term  as  given  by  Eqn.  (25),  however,  does  not  include  the  modification  to  the  pressure  dissipation 
term  as  given  by  Eqn.  (27).  It  can  be  observed  from  Figure  9  that  this  formulation  produces 
fluctuations  near  the  contact  surface  similar  to  what  was  seen  from  the  blended  flux  difference 
procedure.  Additional  calculations  were  done  with  larger  pressure  ratios  that  increase  the  flow  Mach 
number  and  it  was  found  that  the  oscillations  dissipate  when  the  pressure  difference  is  large  enough  to 
generate  a  Mach  number  of  the  order  of  0.1  at  the  contact  interface.  This  provides  conclusive 
evidence  that  this  numerical  instability  is  related  to  acoustic  propagation  at  low  Mach  numbers.  In 
particular,  these  oscillations  appear  to  be  a  fundamental  manifestation  of  inadequate  dissipation  in  the 
pressure  flux  that  arises  for  low  Mach  number  flows  where  the  velocity  arises  from  the  propagating 
pressure  pulse  and  appear  for  both  blended  flux  difference  and  AUSM  family  of  schemes. 
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(a)  Pressure  profile 

Figure  9:  Shock  tube  solutions  for  the  AUSM+up  schemes  with  Mu=1.0 


The  derivation  of  the  additional  pressure  dissipation  in  Eqn.  (27)  was  in  fact  motivated  by  this  test 
case  which  pointed  out  the  inadequacy  of  the  pressure  dissipation  term  in  both  flux  differenced  and 
AUSM  schemes.  The  velocity  dissipation  term  in  the  pressure  flux  for  the  AUSM+wp  formulation 
goes  to  zero  as  the  Mach  number  becomes  low.  However  this  dissipation  may  be  required  for 
reducing  the  pressure  fluctuations  for  low  Mach  number  acoustic  waves  and,  therefore,  the  unsteady 
Mach  number  parameter  was  included  in  Eqn.  (27)  for  the  definition  of  fa.  In  this  so-called 
AUSM+i/  ’p  ’  scheme,  the  addition  of  the  unsteady  Mach  number  scale  in  the  velocity  dissipation  term 
means  that  it  retains  a  substantial  value  for  unsteady  low  Mach  number  flows.  The  effect  of  this 
dissipation  is  illustrated  in  Figure  9  as  indicated  by  the  profiles  given  by  the  green  line.  The  addition 
of  the  velocity  dissipation  allows  for  solutions  that  are  smooth  in  both  the  pressure  and  Mach  number 
profiles.  Furthermore,  the  convergence  is  slightly  improved  and,  therefore,  provides  a  scheme  that 
works  in  an  accurate  and  efficient  manner. 

In  summary,  simulations  for  low  Mach  number  shock  tube  configurations  which  represent  a  pure 
acoustic  wave  propagation  that  generates  a  very  low  Mach  number  flow  indicates  that  both  blended 
flux  difference  and  the  original  AUSM+up/SLAU  schemes  have  some  fundamental  deficiencies. 
While  the  location  of  the  acoustic  front  was  accurately  captured,  the  pressure  and  velocity  in  the 
contact  interface  showed  fluctuations  that  were  a  function  of  the  local  Mach  number.  For  very  weak 
acoustic  waves  where  the  flow  Mach  number  was  0.0001  substantial  oscillations  were  observed  which 
dissipated  as  the  Mach  number  rose  to  0.1  (for  a  larger  initial  pressure  difference).  The  fundamental 
deficiency  in  the  dissipation  was  rectified  in  the  AUSM+u’p’  scheme  by  adding  a  pressure  dissipation 
that  is  enforced  only  for  unsteady  low  Mach  number  flows  and  represents  a  fundamental  advancement 
to  the  AUSM  family  of  schemes.  Unfortunately  for  blended  flux  difference  schemes  there  is  no  clear 
remedy  apart  from  going  to  a  pure  unsteady  preconditioning  procedure  (i.e.  no  blending)  which 
obviously  is  not  acceptable  for  the  broader  class  of  unsteady  low  Mach  number  flows  that  were  tested. 


4.  Concluding  Remarks 

A  generalized  preconditioning  framework  that  is  both  accurate  and  efficient  for  unsteady  low  Mach 
number  flows  is  presented.  It  rectifies  the  deficiencies  of  standard  steady/unsteady  preconditioned 
flux  difference  schemes  for  simulating  unsteady  low  Mach  number  flows  and  has  been  shown  to 
apply  to  a  broad  class  of  problem  encompassing  both  hydrodynamic  and  acoustic  driven  unsteady 
flows.  Two  classes  of  schemes  were  considered:  flux  difference  schemes  with  blended  steady  and 
unsteady  preconditioning  and  AUSM  family  of  flux-splitting  schemes.  For  flux  difference  schemes, 
generalized  “blending”  methodologies  were  developed  (extending  earlier  work  by  Potsdam  et  al.  [1]) 
wherein  “unsteady”  preconditioning  is  used  for  the  pressure  wave  propagation,  while  “steady” 
preconditioning  is  used  for  the  convected  scalars  that  propagate  at  the  fluid  velocity.  Two  well- 
known  AUSM  schemes  were  analyzed  initially;  the  AUSM +up  by  Liou  [5]  and  the  SLAU  scheme  by 
Shima  and  Kitamura  [7]’[8].  The  SLAU  scheme  was  shown  to  be  equivalent  to  the  standard  AUSM+ 
formulation  (i.e.,  without  the  additional  velocity  and  pressure  dissipation  proposed  for  low  Mach 
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number  flows).  A  more  generalized  formulation,  the  AUSM +up’  scheme,  that  provides  a  unified 
framework  for  unsteady  and  steady  flows  using  an  unsteady  Mach  number  parameter  was  developed; 
the  unsteady  Mach  number  parameter  which  can  be  specified  as  a  local  function  of  the  flow  physics  is 
shown  to  provide  the  right  selection  of  dissipation  in  the  mass  flux  term.  Furthermore  a  new  velocity 
dissipation  term  was  developed  in  the  AUSM+t/’/L  formulation  for  acoustic  flows  to  suppress 
spurious  oscillations  at  very  low  Mach  numbers. 

Both  the  blended  flux  difference  schemes  and  the  modified  AUSM+w/?'  schemes  were  tested  over  a 
wide  range  of  unsteady  test  cases  that  encompass  both  hydrodynamic  and  acoustic  unsteadiness.  Both 
schemes  gave  superior  results  with  both  the  pressure  wave  and  velocity  propagation  being  captured 
accurately  independent  of  time  step  and  Strouhal  number  while  providing  adequate  convergence  for 
the  inner  iteration.  However  there  were  two  notable  exceptions:  1)  for  multi-dimensional  flows  such 
as  the  vortex  propagation  problem  the  blended  flux  difference  schemes  distorted  the  pressure  field 
while  the  AUSM +up  ’  scheme  did  not,  and  2)  for  small  pressure  jump  shock  tube  problem,  oscillations 
in  the  contact  surface  were  suppressed  in  the  AUSM+i/  ’p  ’  scheme  by  the  addition  of  a  velocity 
dissipation  term  while  no  remedy  was  available  for  the  blended  flux  difference  procedure.  Therefore, 
the  modified  AUSM +up’  and  AUSM+w ’p’  schemes  are  considered  the  superior  scheme  for  unsteady 
low  Mach  number  flows.  We  note  here  that  the  additional  velocity  dissipation  term  in  AUSM+t/  ’p  ’  is 
needed  only  for  suppressing  the  oscillations  in  the  shock-tube  problem.  In  all  other  cases,  the 
modified  AUSM+w/U  scheme  proved  to  be  the  best  choice.  Resolving  this  discrepancy  will  be 
addressed  in  future  work.  Additionally,  improvement  of  the  inner  iteration  convergence  for  the 
AUSM+wp'  scheme  also  deserves  further  scrutiny.  Here,  the  convergence  may  currently  be  hindered 
by  the  inconsistency  between  the  LHS  (implicit  matrix)  and  RHS.  Finally,  future  work  will  also  focus 
on  extending  the  methodology  to  frameworks  that  have  both  rotational  and  inertial  frames  as  well  as 
on  developing  generalizing  the  unsteady  preconditioning  parameter  to  consider  local  and  global 
Strouhal  numbers. 
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